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Abstract 

We study the phenomenon of lack of reversibility in molecular dy- 
namics algorithms for the case of Wilson's lattice QCD. We demon- 
strate that the classical equations of motion that are employed in 
these algorithms are chaotic in nature. The leading Liapunov ex- 
ponent is determined in a range of coupling parameters. We give 
a quantitative estimate of the consequences of the breakdown of 
reversibility due to round-off errors. 
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1 Introduction 



Ever since Wilson's formulation of QCD for a Euclidean lattice, nu- 
merical simulations have played an important role in non-perturbative 
investigations of physics concerning the strong interaction. When per- 
forming numerical simulations of lattice QCD, especially with dynamical 
fermions, one usually relies on molecular dynamics kind of algorithms 
such as the Hybrid Monte Carlo (HMC) [|I| - the most popular one - 
or the Kramers equation algorithm g, g]. 

A major ingredient of these algorithms is the -numerical- integration 
of a set of non-linear differential equations, the classical equations of 
motion. The non-linearity of these equations may give rise to the sus- 
picion that the equations of motion may behave like those of a chaotic 
classical dynamical system. Indeed, it was first shown in and later 
confirmed in [g, ^ that the equations of motion used in the HMC 
algorithm, when applied to QCD, are chaotic in nature. This amusing 
but, at first sight, academic observation has an important consequence 
for any practical application of such kind of algorithms. 

The exactness of the molecular dynamics kind of algorithms is guar- 
anteed by the -sufficient- detailed balance condition. For this condition 
to hold, the equations of motion for the system, which govern the clas- 
sical motion, are to be reversible. However, as noticed some time ago 
[[7j], this reversibility condition is violated due to round-off errors occur- 
ring in the numerical integration of the equations of motion. The fact 
that the equations of motion are chaotic in nature now means noth- 
ing else but that the round-off errors get magnified exponentially along 
the integration of the equations of motion with a positive Liapunov 
exponent v. This by now well-established phenomenon p, 0, 0, §, || 
therefore leads to rounding-error effects much larger than might naively 
be expected. Because of these violations of the reversibility condition, 
the algorithms are no longer exact from a principle point of view. Of 
course, the question remains at what quantitative level the reversibility 
violation manifests itself in physical observables. 

In this paper, we provide a detailed study of the Liapunov exponents 
for the case of QCD and determine them in a range of values for the 
couplings of the theory. A short account of our results has appeared 
in [^. The dependence of the leading Liapunov exponent on the cou- 
pling constants allows us to estimate where, in the coupling parameter 
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space, we may expect a stronger or weaker exponential amplification 
of rounding-error effects. We then try to estimate how much the re- 
versibility violation affects physical observables such as the plaquette, 
the Polyakov line, plaquette correlation functions and the scalar density. 
Finally we emphasize the special role of the acceptance Hamiltonian 
with respect to the rounding-error effects. 

Although we will discuss these effects here for the case of lattice 
QCD, we think that similar problems may also arise in different fields 
where molecular dynamics algorithms are used, e.g. in fluid dynamics. 
There again it might be that the equations of motion used resemble 
those of a chaotic, classical, dynamical system with a positive Liapunov 
exponent. In some cases, as in part of the QCD simulations presently 
done, computers with 32-bit arithmetic are used. For these case, espe- 
cially, it is desirable to know whether round ing-errors have a noticeable 
and significant influence on the values of physical observables. In addi- 
tion, to speed up the numerical simulations, even on 64-bit arithmetic 
machines the time-consuming part of the computational work is some- 
times performed in 32-bit arithmetic, for which case one would again 
be interested in the rounding-error effects. 

2 Molecular dynannics equations for lattice QCD 

Let us briefly sketch the HMC algorithm that is used for simulations of 
lattice QCD. Consider first for simplicity a quantum mechanical system 
quantized by Feynman's path integral prescription and specified by its 
action S{q), where q denotes a path of the quantum mechanical particle. 
In order to evaluate Feynman's path integral numerically, one sets the 
system on a Euclidean discrete time lattice. The task is then to gen- 
erate classical paths q that are distributed according to the Boltzmann 
distribution e"'^. 

One way to proceed is to use a method that resembles those used 
in molecular dynamics simulations. To this end one introduces an ad- 
ditional, fictitious so-called Monte Carlo time t and corresponding mo- 
menta p. The HMC algorithm then works as follows. 

One starts from some given path qi on the discrete time lattice and 
generates initial momenta ^ from a Gaussian distribution of unit vari- 
ance and zero mean. From this initial set {qt,pi) one computes an initial 
classical Hamiltonian H{pi,qi), where h= \p^ + s{q). The fictitious time 
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evolution of the momenta and the paths is now given by the following 
set of coupled first-order differential equations: 

p = — T- , q^p ■ (1) 

dq 

The time derivatives in eqs. (|I]) are to be understood with respect to the 
fictitious Monte Carlo time. The numerical integration of the equations 
given in eqs. (p is performed by using a discretized form of eqs. (|T]). In 
practice one uses a leap-frog integration scheme, using Nmd integration 
steps of size 6t in order to integrate from fictitious time to some value 
r. Such an integration procedure is called a trajectory of length t. 

In order to generate the desired Boltzmann distribution and to ac- 
count for the discretization errors of the integration procedure, the new 
momenta pf and coordinates g/ obtained after the integration of the 
equations of motion are only accepted with a certain probability P, 
p = min(i,exp{-Ai/}) with AH = H{pf,qf ) - H{pi,q,). The Hamiltonians 
that enter this Metropolis step are called the acceptance Hamiltonians. 
The desired distribution is now obtained by repeating this procedure 
a sufficiently large number of times. The above prescription ensures 
the detailed balance condition, and the resulting algorithm is hence ex- 
act in the sense that it indeed converges to the anticipated Boltzmann 
distribution. 

The heart of the algorithm sketched above is the integration of the 
classical equations of motion eqs. (0). This part is reminiscent of tech- 
niques used in molecular dynamics algorithms. Quite often the equa- 
tions of motion derived from the problem under consideration are non- 
linear and may as such correspond to those characterizing a chaotical 
dynamical system. 

Let us now come to our example of lattice QCD, which is a relativis- 
tic field theory defined in 4-dimensional Euclidean space-time. To be 
specific, we will focus our discussion on the equations that arise in the 
simulations of QCD in Wilson's formulation with two flavours of quarks 
with degenerate masses. 

Lattice QCD is established on a Euclidean space-time lattice of size 

X T. With lattice spacing set to unity, a point x on the lattice has 
integer coordinates x = {t,xi,x2,x3), which are in the range < i < r;0 < 
x^ < N. A gauge field U^^f, e SU{3) is assigned to the link pointing from 
point X to point + where ^ = 0,1,2,3 designates the four forward 
directions in space-time. In this study, periodic boundary conditions 
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have been taken for the gauge fields and for the quark fields in all four 
directions. The full partition function for Wilson QCD is given by 

Z = JvUexp{-Sg)T)ct{Q^) = JvUV(j)^V(f>exp {-Sg-(t)^Q-^(l)) . (2) 

The term Sg in the exponential is the pure gauge action and is given by 

5, = -^5]Tr(L/p + 4) . (3) 

p 

The symbol Up represents the usual plaquette term on the lattice. 
The determinant factor Det{Q^) represents the contribution of internal 
fermion loops to the theory. The matrix Q that appears in the determi- 
nant is a Hermitian sparse matrix defined by: 



Q{U)x,y = Co75 



5x..y -l^^i^ - lt.)U^,t.Sx+f,.y + (1 + lt^)Ul^,,.j,Sx-^,,y 



(4) 



with K the so-called hopping parameter and cq = i/(i + 8k). The aim of 
the simulation is to generate configurations according to the probability 
distribution exp{~Seff) = exp{~Sg - (j)^Q^'^(t)) using Monte Carlo methods. 
The Hamiltonian for the simulations of lattice QCD is given by 



X.fl 



where i/^;,^ is the momentum conjugate to the gauge field u^^f, and takes 
values in su(3), the Lie algebra of S';7(3). 

Since the bosonic part Sb = (t)^Q~'^(l) is quadratic in the ^ fields, these 
are generated at the beginning of each molecular dynamics trajectory via 
(j) = QR, where i? is a random spinor field which is Gaussian distributed. 
The kinetic term in eq. (§) is also generated from Gaussian noise at the 
beginning of the update. Then the gauge fields and their corresponding 
momenta are updated according to the equations of motion: 

Ux.fi = iHx,iiUx,i_i , iHx,i_i = \Ux,i_iFx.fi\T.A. , (6) 

where the symbol [■■■]t.a. stands for taking the traceless anti-Hermitian 
part of the matrix and the quantity Ux,^Fx^f, is the total force as- 
sociated with the link Ux,,,. The dot on a field variable represents the 
derivative with respect to the Monte Carlo time t. The quantity f^.^^ is 
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nothing but the "coefficient" in the change of the effective action when 
an infinitesimal change of the gauge link su^^f, is applied, i.e. 

SS^ff ^Y.^t{F,,^SU,.^ + FljUl^) . (7) 

XfJ, 

One can easily check that the time evolution described by eqs. (|5|) 
conserve the Hamiltonian (^. Equations (|5D define a Hamilton flow in 
a phase-space manifold that is a direct product of m^T factors of su{3) 

and su(3). 

Equations (6) establish the molecular dynamics part of the HMC 
algorithm and it is in this part where the problems with rounding- 
errors appear. Note that one considers here a purely classical set of 
autonomous first-order differential equations. 



3 Liapunov exponents 

In this section, we briefly describe the concept of the Liapunov charac- 
teristic exponents [|ID|, [TI|, [T^. Liapunov exponents serve as an impor- 
tant quantitative measure for the degree of stochasticity of a dynam- 



ical system. For an introduction to this topic, see [[T^ and references 
therein. 

Let us consider a time evolution g(r), p(r) , described by a set of first- 
order autonomous differential equations like eqs. To each point 
in the phase-space manifold {q,p) one can attribute a local Liapunov 
exponent. These exponents describe the mean exponential rate of di- 
vergence of the distance between the trajectories of two nearby points in 
phase-space. Given two points in the phase-space manifold, which are 
close to each other at r = o, one can follow the time evolutions of the 
trajectories originating from these two points. At any given instance of 
the time r, one can construct the vector pointing from a point on one 
of the trajectories to the corresponding one on the other trajectory. If 
the distance between the two initial points in phase-space becomes in- 
finitesimally small, the above-mentioned vector belongs to the tangent 
space of the phase-space manifold. 

The Liapunov exponents can now be determined by studying the 
time evolution of these tangent vectors. In order to be specific, we 
will directly discuss the example of the equations of motion used in the 
HMC algorithm as applied for lattice QCD. We start denoting by dH^.j, 
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and dX^^f, = -iU-ldU^^f, the tangent vectors of H^^f, and 17^^,,, respectively. 
The time evolution of these vectors is described by 



dXx.p, = U^^^dHx^fiUx.^j. , idHx.^j, — d\Ux.^j.Fx,n]T.A. , (8) 

which follow from eqs. One can study this set of equations for 
{dHx^f,,dXx,^) together with the original eqs. (§). Note that eqs. (|8|) are 
linear in dH^^f, and dx^,,,, both of which are elements of su(3). Let us 
now assume that, initially at r = o, {dHx,f,,dXx,f,) were given some value 
in the tangent space of some point in the phase-space manifold. We 
can introduce the norm in the tangent space as: 

D'{r) = Tr(di?L(^) + dX^J^)) . (9) 

X,fJ, 

As time evolves, this norm will change and the Liapunov exponent can 
be defined as: 

1^ = lim lim — log ^[ } . (10) 

In fact, due to its linear nature in {dHx^^,dXx^^), eqs. (§) can be written 
as 





J[Hx,,,Ux,,]-[ "^''^ , (11) 



and the Liapunov exponents are just the time averages of eigenvalues 
of the matrix J along a trajectory. These exponents can be ordered 
according to the magnitude of their real part. The one with the largest 
real part is called the leading Liapunov exponent and will be studied 
in more detail below. It is known []ri| that Liapunov exponents do not 
depend on the choice of metric in phase-space. 

The numerical calculation of Liapunov exponents of a given flow can 
be done straightforwardly [|I2|. We now go to the discretized versions]] 
of eqs. (§) and eqs. (§) and integrate them simultaneously in time using 
a leap-frog integration scheme. As a safety measure, in order to avoid 
severe rounding-errors, the tangent vector should be renormalized with 
respect to its initial norm after each integration step. The procedure 
goes as follows: Starting at r = o, we have some initial value of d'^{0) = 
Y.x.,,^^{dHl ^ + dxl ^) for a given initial tangent vector {dHx,f,,dXx,f,). We 

*ln the discretized version of numerical integration, we replace the differentials by the 
corresponding differences. 
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then integrate one step in time with step size St using the leap-frog 
scheme. Now, we evaluate the norm d{St) of the new tangent vector and 
store this information for the computation of the Liapunov exponents. 
Next, we rescale the new tangent vector in such a way that its norm 
is still equal to D{0), i.e. (dH,,^^{ST),dX^,^{ST)) = {dH.,,^{5T),dX.,,^{5T))/D{5T). 
We then proceed to the next integration step, starting from this already 
normalized tangent vector. Because of the linear nature of eqs. this 
rescaling is legitimate since the Liapunov exponents only depend on the 
ratio of the norms. We then proceed to integrate over n such steps. It 
can be shown [0 that the average 



Vn. 



1 DikSr) , , 



nSr ^ " D(0) 



is approaching the leading Liapunov exponent when n ^ oo. In addition, 
its value is independent of the value of 5t, as long as the step size 5t is 
not too large. 

In contrast to the truly classical Hamilton flows, for which the phase 
space is usually partitioned into regular and stochastic regions, in the 
case of the HMC algorithm considered here, the exponent we find is a 
phase space average of local exponents weighted with the appropriate 
Boltzmann factor. 

Another remark is that the concept of Liapunov exponents could 
be generalized to higher-dimensional objects as well; not surprisingly, 
the leading exponent for a p-dimensional volume formed by p linear 
independent tangent vectors is the sum of the p leading exponents of 
the tangent vector. We point this out because this observation can 
serve as a tool for calculating the subleading Liapunov exponents if one 
wishes to. 



4 Liapunov exponent for simulations of QCD 

We have used the method described in the previous section to study 
the Liapunov exponents in simulations of QCD. A determination of the 
leading Liapunov exponent for QCD simulations with gauge group SU{2) 
was reported in [^. Similar studies were done for gauge group 5C/(3) in 
[0, ^. Here, we would like to extend these studies by computing the 
Liapunov exponents in a range of coupling parameter values. 
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When numerically integrating the equations of motion derived from 
the lattice QCD action, we have chosen a step size of 5t = 0.02, which 
results in an acceptance of almost 100%. We have integrated the equa- 
tions of motion with 200 steps, which amounts to a trajectory length of 
4.0. The calculations are mainly done on a 4'' lattice but we also used 
8" and 16^ lattices. We checked that our results do not depend on the 
value of the step size for 5t < 0.02. The initial values of {dHx^i,,dXx^f,) are 
generated from a Gaussian distribution. Then the norm of the tangent 
vector is recorded for each time step, from which we could define the 
"instantaneous" Liapunov exponent at the fc-th step as: 

Note that, as discussed above, in the actual simulation l»(t) is always 
normalized to the norm of the starting vector in order to avoid potential 
rounding-errors. In order to monitor the change of diJ^;,^ with t and dx^^f, 
separately, we have also recorded the values of ^'^(r) and Dx{t), which 
are defined as 

DUr) = 5]Tr(di/J^^(r)) , D^ir) = ^ Tr(dX,2_^(r)) . (14) 

XfJ, XfJ, 

We also define the "instantaneous" -or effective- exponent for dH^^^, as: 

^H(r) = ^\og^^ (15) 
St Dh(0) 

and similarly i^xir) for dx^^,^. Again we remark that Dnir) receives a 
proper normalization in the simulation. 

We first describe our results for the Liapunov exponents obtained 
in pure SU{3) gauge theory on 4"^ lattices. We have studied the theory 
for various values of the gauge coupling (3 ranging from /? = 0.5 up to 
^ = 30. In fig. |I|, we show the behaviour of the instantaneous Liapunov 
exponents as a function of the trajectory length for various values of 
f3. From these figures, it is clearly seen that the exponents show an 
oscillatory behaviour, especially at the beginning of a trajectory. Then 
the amplitude of the oscillation dies out as time evolves. Finally, all 
three exponents u, vh and vx approach the same stable constant value. 
The leading Liapunov exponent can be extracted by filtering out the 
zero frequency part in v{t) towards the end of the trajectory. In our 
determination of the leading exponent, we fit the tail of the function 
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Figure 1: The instantaneous Liapunov exponents, as defined in eq. ( p!3| ) and eq. (|T^), 
are plotted as a function of the IVlonte Carlo time t along a trajectory for various 
values of (3. The dotted lines represent vh{t) while the dashed lines represent vx{t). 
The v{t)'s are represented by solid lines. 



v{t) to a constant. We also note that as the value of I3 is increased, so 
is the frequency and magnitude of the oscillation. This is a reflection of 
the common expectation that for large values of (3, the system is getting 
closer to a Gaussian model. Despite rather strong oscillations for v{t) 
at large values of/?, the integral r{t) of the function v{t), defined as 

fc 

R{t) = ^ 5Tv{j8T), T = kSr, (16) 

which reflects the change in the norm relative to the starting point, 
is basically increasing linearly with only little oscillations, as shown in 

fig- @- 

Having extracted the leading Liapunov exponents [| for several values 
of /3, we plot them as a function of (3 in fig. ^ (a). We see from this 
figure that the Liapunov exponents show a significant /3-dependence. 

t Similar studies have also been done by the authors of ref . [|j . Complete consistent values of the 
exponent have been obtained, except for very large (3 values. 
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12 3 4 



Figure 2: The integral of the function i^{t), eq. (|16|), which represents the change 
in the norm of the tangent vector as compared to its original value, is plotted as a 
function of the Monte Carlo time of a trajectory for the same values of (3 as in fig. |^. 



From /? = 5.0 to /? lO.O, we have a leading Liapunov exponent that 
ranges from 0.85 to 0.35. 

In fig. 0(a), we plot also results for the Liapunov exponents for an 
8" (filled triangles) and for a 16^ lattice (stars). We can see that the 
Liapunov exponents for the different size lattices fall essentially on one 
common curve. This behaviour is not in accord with the anticipated 
property of the Liapunov exponent to be inversely proportional to the 
correlation length, as proposed in [Q. The data seem to indicate that 
we are really studying here classical equations of motion derived from 
the 4-dimensional classical Hamiltonian given in eq. (||). 

We have also performed similar investigations of the Liapunov ex- 
ponents for full QCD with dynamical fermions. In dynamical QCD 
simulations, the parameter k appears and the value of the Liapunov 
exponent, in principle, will also depend on this parameter. For k val- 
ues that we have studied on 4^ lattices, the leading Liapunov exponent 
shows, however, only a rather weak dependence on n. 

We have performed a rough scan in the parameter ranges 5.6 < /? < 8.0 
and 0.13 < K < 0.16. The Liapunov exponents are measured in the same 
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Figure 3: The average leading Liapunov exponents for pure SU{'i) gauge theory (a) 
and full QCD (b) as a function of the gauge coupling /3. In (a), the open and filled 
triangles correspond to 4* and S** lattices, respectively. The stars denote our results 
for a IG** lattice. In (b), exponents for different values of n are represented by different 
symbols. The triangles correspond to k = 0.13 and the hexagons to k = 0.16. 



way as in the pure gauge case. The instantaneous exponents also show 
similar behaviours. In fig. 0(b), we show the final result of the exponents 
as a function of /? and k. Different n values for the same /3 give almost 
the same value for the Liapunov exponents. For each value of (3, we 
have taken four different values of k, namely o.i3, o.i4, o.i5 and 0.16. 
Due to their weak dependence on k, we only show the exponents for 
two values of k, i.e. k = 0.13 (triangles) and n = 0.I6 (hexagons). The 
errors of the points are about the size of the symbols. In the parameter 
range that we have studied, the Liapunov exponent is roughly 0.4-0.5. 

5 Consequences of irreversibility 

The equations of motion in eqs. (^) are to be reversible to ensure the 
exactness of the HMC algorithm used for the simulations of lattice 
QCD. However, rounding-errors in the numerical integration process 
of this algorithm imply violations of the reversibility condition. Thus, 
the exactness of the algorithm can no longer be guaranteed. Although 
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this brings us into an unsatisfactory situation from a principle point of 
view, it is necessary to investigate what are the consequences of the 
reversibility violations. 

For a molecular dynamics trajectory there are two causes for an in- 
crease of rounding-errors when the trajectory length r is increased for 
fixed integration step size St. If we assume that the rounding-errors are 
Gaussian distributed, we expect a random walk behaviour and hence 
that the rounding-errors grow like ^/T. Second, as demonstrated above 
and also found in [§, 0, ||, ||, since the equations of motion behave like 
those of a classical chaotic system, the rounding-errors get amplified 
exponentially, like exp(i/T), with a positive Liapunov exponent u. 

Both effects lead to a certain amount of reversibility violation and 
in this section we try to quantify the effect of this reversibility violation 
on physical observables on the one hand and on the Hamiltonian that 
enters the Metropolis decision in the HMC algorithm on the other hand. 

5.1 Observables 

To study the effects of the above discussed reversibility violation on 
observables, we have run two HMC simulation programs in parallel. 
One of them is with 32-bit arithmetic but with 64-bit arithmetic scalar 
products and summations over the lattice, which mimics the typical 
situation in simulations on a 32-bit machine. The other one runs with 
complete 64-bit arithmetic and serves as a reference point for an "exact" 
program. Let us remark that for this section we have chosen SU{2) as 
the gauge group. 

We then proceeded as follows: one generates a configuration, say 
with the 64-bit arithmetic program version. On this configuration one 
starts the 32-bit and the 64-bit program versions and run them for 
a given number of molecular dynamics steps with a fixed step size, 
i.e. up to some trajectory length r. At each step of the trajectory, 
one measures quantities such as the plaquette P, the Polyakov line L, 
plaquette correlation functions or, in the case of dynamical fermions, 
the scalar density. Let us denote such a quantity with O for the 64-bit 
version of the program and with d the corresponding one for the 32-bit 
arithmetic version. As an example, let us give the measurement of the 
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Polyakov line: 

1 ^"^ 

t=o 

i.e. the trace of the ordered products of gauge field variables. The 
index f denotes the 3-dimensional space coordinate. We then measure 
at each molecular dynamics step 

\\L-Lr^^j2(Ls^Lsr. (18) 

X 

For other observables the difference \\o - d\\ can be computed in an 
obviously generalized analogous way. The whole procedure is repeated 
on a number of configurations in order to obtain an error estimate. 
From the above discussion we expect that 

\\0-d\\ =co + ciV^ + C2P/" . (19) 

The result for the Polyakov line on a 6^8 lattice is shown in fig. 0, 
together with a fit to the data according to eq. ([T^ ) for both the pure 
gauge theory (open squares) and dynamical fermion simulations (filled 
squares). Obviously, the fit formula eq. ([T^) provides an excellent de- 
scription of the data. We checked that also the other above-mentioned 
observables can be fitted by eq. (|I^). We remark that for the scalar 
density the exponential growth of round ing-errors sets in for only large 
trajectory lengths, t « 8. We conclude that indeed both sources of 
rounding-errors, the random walk and the chaotic behaviour, seem to 
be present in our simulations in practice and show up in physical ob- 
servables. 

We find the same rounding-error behaviour of eq. ( \[% when we first 
sum up the observables over the lattice, constructing Og = Y^xOx, and 
then build the g/o&a/ difference \Og-dg\ at each molecular dynamics step. 
The prominent exception is the global difference of the plaquette, \Pg- 
Pg\, which shows an oscillatory behaviour similar to the ones shown in 
fig. 1. The reason for these oscillations is that in the HMC algorithm the 
total Hamiltonian is, up to finite step size errors, conserved. Therefore, 
the tangent vector that corresponds to the plaquette has to compensate 
the oscillations of the tangent vector that corresponds to the momenta 
studied in section 4. 

For a lattice size of 6^8 we see from fig. @|that the difference |li-L|l, 
which indicates the deviation from the true value, is very small and hard 
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Figure 4: The difference \\L — L\\, eq. (18), for the Polyakov line L as a function of 
the trajectory length r. Open squares are for the pure gauge theory at (3 — 2.12 and 
filled squares are from full dynamical fermion simulations at = 2.12, k = 0.15. The 
small inset amplifies the region of small r. The solid line is a fit according to eq. (19). 



to detect in a simulation. By true value we here mean -by definition- 
the value of the observables from the 64-bit arithmetic program version. 
It is only for very long trajectories that one would risk the danger to 
produce significant differences. 

One may expect that the coefficients c, in eq. ([T^) depend on the 
lattice size. Our original plan was therefore to extract the volume de- 
pendence of these coefficients and to extrapolate to a lattice size, where 
rounding-error effects and in particular the consequences of irreversibil- 
ity will become dangerous for trajectories of length i. Our findings for 
the coefficients for on various lattice sizes is given in table |I|. 

To our surprise, the fit parameters show basically no lattice size depen- 
dence at all. A similar finding was obtained for the plaquette variable 
and for plaquette correlation functions at various distances. 
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Table 1: Lattice dependence of the fit parameters in eq. (|l9|) for the Polyakov line 
at /3 = 2.12 for the quenched case and (3 = 2.12, k ~ 0.15 for the unquenched 
simulations. 



Lattice 




ci • 10^ 


C2 • 106 


V 


Co • 106 


4^8 


quenched 


0.1180(10) 


0.108(2) 


0.630(6) 


-0.027(2) 


6^8 


quenched 


0.1177(6) 


0.104(1) 


0.654(3) 


-0.023(1) 


84 


quenched 


0.1177(6) 


0.103(1) 


0.659(3) 


-0.022(1) 


4^8 


dynamical 


0.1111(14) 


0.143(4) 


0.557(8) 


-0.062(4) 


6^8 


dynamical 


0.1066(21) 


0.143(5) 


0.608(10) 


-0.060(5) 


84 


dynamical 


0.1115(17) 


0.134(4) 


0.623(8) 


-0.052(4) 



For the HMC algorithm applied for the pure gauge theory one may 
give a simple reasoning: since in the HMC algorithm the gauge links 
only feel their nearest neighbours, their update is essentially local. But, 
for dynamical fermions, in the equations of motion the information of 
the inverse fermion matrix is used, which is non-local. The results of ta- 
ble |I] are therefore somewhat counter-intuitive. This indicates that the 
difference between the 32-bit and 64-bit versions in the solution vectors 
obtained from the conjugate gradient (CG) method do not increase 
significantly with the volume. We checked this explanation explicitly 
by studying the stability of the CG inversion method against rounding- 
errors. In addition it might be that our studies are performed in a situa- 
tion where the propagators vanish fast with growing distance. Therefore 
information would only be transported over a few lattice spacings. 

5.2 Hamiltonian 

Would the HMC algorithm only contain the molecular dynamics part, 
the results of the previous section would indicate that rounding-error 
effects show up in physical observables only for a very large number 
of integration steps. However, the exactness of the HMC algorithm 
requires a Metropolis reject/accept step for which the difference of the 
initial and final values of the Hamiltonians are taken. An important 
observation is that for this difference an absolute precision is required. 
The computations of the acceptance Hamiltonians therefore play a spe- 
cial role and an investigation of rounding-error effects on them is most 
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crucial. 

We therefore measured in full QCD the difference between the value 
of the Hamiltonian H from the 32-bit arithmetic program version and 
the one from the 64-bit arithmetic program version H. Let us define 
the ratio 

_ (l^(Ag)l) 

" (|Ai7|) ' ^^^^ 

where {\6{/^H)\) = (|AH-AiJ|) and ah is the difference between the initial 
Hamiltonian at the beginning and the final Hamiltonian evaluated at the 
end of a trajectory. One should expect a dependence on the lattice size 
for Rh. Keeping the acceptance rate constant means that (|Ai/|) should 
be constant. However, the values for H themselves increase substantially 
with the lattice size since they are constructed by a sum over all lattice 
points, see eq. (^. As a consequence, the rounding-errors become more 
and more significant relative to (|AiJ|). 

Indeed, averaging over 100 trajectories, we find Rh ~ o.08%,o.i8% 
and 0.70% for 8"*, 12^ and 16^ lattices, respectively. These numbers are 
obtained by performing only 10 molecular dynamics steps with step size 
5t = 0.02, a number that is certainly too small for simulations on a 16'' 
lattice. Still, extrapolating Rh linearly in N"^ to a 32'' lattice we find that 
(|(5(Ai?)|> could become about 5% of (|Ai?|). We think that this number is 
a conservative estimate of the rounding-error effects on the acceptance 
Hamiltonian on such a large lattice. We feel that performing simulations 
on lattices of comparable sizes on 32-bit arithmetic computers could not 
be safe. The same holds, of course, on 64-bit arithmetic machines when 
32-bit arithmetics is used for the numerical computations. 



6 Conclusions 



In this paper, we have investigated several questions related to the 
reversibility problem of the molecular dynamics kind of algorithms as 
used for the simulation of lattice QCD. The purely classical equations of 
motion used in these algorithms behave like those of a chaotic dynamical 
system. We have determined the leading Liapunov exponents v for both 
pure SU{i) gauge theory and full QCD for various bare parameters. We 
have found a weak dependence of the leading Liapunov exponent on 
the hopping parameter «, whereas the /3-dependence is more significant. 
We argued that, due to the chaotic nature of the equations of motion. 
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the rounding-error that occurs in the integration of these equations is 
magnified exponentially with increasing trajectory length, when keeping 
the step size constant. 

We suggested that the growth of the round ing-errors as a function 
of the trajectory length r should consist of two parts: the first is a 
random walk behaviour, leading to a ^ part. The second is the ex- 
ponential amplification, leading to an exp(z/T) contribution. We verified 
that this ansatz is indeed followed by observables such as the plaquette, 
the Polyakov line and plaquette correlation functions. Considering the 
above observables, the rounding-errors, however, do not increase with 
the lattice size, at least for the parameter values we considered here. 
The situation is different for the value of the difference AH of the initial 
and final Hamiltonians, computed at the beginning and the end of a 
molecular dynamics trajectory, which enters the Metropolis decision in 
the HMC algorithm. This quantity plays a special role in the molecular 
dynamics kind of algorithms, because it is the absolute precision with 
which these Hamiltonians have to be calculated. We find that, with 
growing lattice size, the rounding-errors for AH increase. Extrapolating 
the rounding-error effects as found on 8*, 12* and I6* lattices to a lattice 
of size 32^, gives a rounding-error for AH that can reach several per cent. 
We regard this estimate as rather conservative and therefore conclude 
that simulations on lattices of size 32'' using 32-bit arithmetics, either 
by hardware or by software, could be dangerous. 

A real test of the effects of rounding-errors on observables would be 
to run a 32-bit arithmetic program version against a 64-bit arithmetic 
version for a long time and see, whether one finds differences in some 
observables. Here we could only perform approximations of this test. 
We can therefore not give a definite answer of whether rounding-errors 
will lead to problems for the molecular dynamics kind of algorithms. 
However, our data indicate that it can be dangerous to use these algo- 
rithms for simulating large lattices employing 32-bit arithmetics. 
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